Optimal extraction of data contained within a signal requires the removal of defects such as noise and instrumental limitations. A key area in which optimized extraction and reconstruction of data is sought is the field of image enhancement. Even when instruments can be made virtually noise-free, instrumental signatures related to finite spatial, spectral, or temporal resolution remain. At this point, image reconstruction is required to remove the instrumental signatures. Applications of image enhancement, and the sources of noise and other factors that can negatively impact data extraction, cover a wide range including astronomical observation and planetary exploration, where sources can be faint and atmospheric interference introduces noise and distortion, military and security surveillance, where light can be low and rapid movement of targets result in poor contrast and blur, medical imaging, which often suffers from lack of clarity, and video images, where transmission and instrument limitations, and the need for real time response, can negatively impact image sharpness and detail.
Digital image processing has been developed to provide high-quality, robust reconstructions of blurry and noisy data collected by a variety of sensors. The field exists because it is impossible to build imaging instruments that produce consistently sharp pictures uncorrupted by measurement noise. Nonetheless, it is possible to mathematically reconstruct the underlying image from the nonideal data obtained from real-world instruments so that information present but hidden in the data can be extracted with less blur and noise. Many such methods utilize a process in which a predictive model is constructed and compared to the data to evaluate the accuracy of the model's fit to the data.
Throughout this written description, “data” refers to any measured quantity, from which an unknown “image” is estimated through the process of image reconstruction. The term image denotes either the estimated solution or the true underlying image that gives rise to the observed data. The discussion usually makes clear which context applies; in cases of possible ambiguity “image model” is used to denote the estimated solution. Note that the data and the image need not be similar and may even have different dimensionality, e.g., tomographic reconstructions seek to determine a 3D image from projected 2D data. An alternative term to image is “object”, which conveys the idea that the model can be more general than an image. The two terms are used synonymously hereafter.
Statisticians have long sought to restrict the parameters used to fit data in order to improve the interpretation and predictive accuracy of the fit. The standard techniques are subset selection, in which some of the parameters are found to be unimportant and are excluded from the fit (e.g., Miller 2002), and ridge regression, in which the values of the parameters are constrained by adding a regularization term to the merit function used in the fit (e.g., Tikhonov 1963). Tibshirani (1996) combined the two methods in a technique known as Least Absolute Shrinkage and Selection Operator (LASSO).
The need to limit the number of parameters in a fit is imperative for underdetermined or poorly determined problems, in which the number of parameters is greater than or comparable to the number of data points. As already emphasized by the originators of the statistical methods of minimum least squares (Gauss 1809) and maximum likelihood (Fisher 1912, 1922), these methods are only valid in the asymptotic limit in which the number of data points greatly exceeds the number of fitted parameters. Away from this asymptotic limit, noise is fitted as signal, and the fit loses its interpretative and predictive power.
Puetter er al (2005) review numerous reconstruction algorithms in current use, including iterative image reconstruction methods, which iteratively fit image models to the data. Many prior art iterative schemes that are designed to converge to the maximum-likelihood solution converge slowly, even when they are terminated early to avoid overfitting. Convergence can be achieved more quickly by using the Hessian matrix of second-order partial derivatives of the merit function with respect to the variable (Hesse 1876). Unfortunately, even using this approach, the Hessian matrix is too large to be computed when dealing with the large-scale problems that are frequently encountered in image reconstruction, e.g., where the image contains a large number of pixels with significant emission. In such cases, the matrix elements can number in the trillions. Matrices of this size simply cannot be processed by present-day computers, or even stored in memory.
The desire to restrict the number of parameters in underdetermined or poorly determined problems—sparsity, as it is called today—is based on the more abstract concept of minimal complexity (Solomonoff 1964; Kolmogorov 1965; Chaitin 1966), which goes back to the medieval work of William of Ockham, who advocated parsimony of postulates. Put simply, other things being equal, a simpler explanation is better than a more complex one.
Finding a sparse solution, e.g., by adding an l0 norm regularization term to the merit function, is an N-P hard problem, in which the computational effort increases more rapidly than any polynomial in the number of parameters. This has led to the replacement of the l0 norm with an l1 norm (Chen, Donoho & Saunders 1999), so that the optimization of the fit with respect to the parameters becomes a solvable convex problem. Candès, Romberg & Tao (2004) went a step further and showed how randomly to reduce the amount of data required for the fit under a condition of incoherence, a technique known as compressed sensing. Donoho (2006) showed that, under a similar condition of incoherence among the basis functions used in the parameterization, the minimal l1 norm solution is also the sparsest solution.
The disadvantages of the l1 norm methods are twofold. First, many problems of interest simply do not satisfy the incoherence condition and are unsuitable for l1 norm methods. Second, even when they do satisfy the incoherence condition, large-scale problems require excessive computational effort. Thus, although they are convex problems and solvable in principle, they cannot be applied in practice to present-day problems with millions or more parameters, a problem that also plagues traditional statistical methods. Donoho et al (2006) discuss how to apply randomness more efficiently to such large-scale problems without the use of an l1 norm.
The pixon method is an efficient technique to obtain minimally complex solutions of pixel-based data—including large-scale problems—without demanding strict sparsity and without the incoherence condition required by methods utilizing the l1 norm (See, e.g., Piña & Puetter 1993; Puetter & Yahil 1999; Puetter, Gosnell & Yahil 2005, and U.S. Pat. Nos. 5,912,993, 6,353,688, 6,490,374, 6,895,125, 6,993,204, 7,863,574, 7,928,727, 8,014,580, 8,026,846, 8,058,601, 8,058,625, 8,086,011, 8,090,179, 8,094,898, 8,103,487, the disclosures of which are incorporated herein by reference). The pixon method is therefore useful for large-scale underdetermined or poorly determined inverse problems such as image reconstruction or spectral analysis. Minimum complexity is achieved by adaptively smoothing at each pixel location by the widest kernel among a library of kernels, such that smoothing with this and all narrower kernels provides an adequate fit to the data footprint of the pixel under consideration. The map specifying which kernel to use at each pixel is called the pixon map.
In its current form, a pixon reconstruction consists of three steps. First, it reconstructs a “pseudoimage” without any pixon constraints. Second, this pseudoimage is used to determine the pixon map. Third, the final image is obtained by a constrained reconstruction guided by the pixon map. Steps two and three can be repeated a number of times, but this is usually not necessary in practice, provided that a reasonable pseudoimage is obtained in the first step. See, Puetter & Yahil (1999), Puetter et al (2005) and U.S. Pat. Nos. 5,912,993, 6,353,688, 6,490,374, 6,895,125, 6,993,204, 7,863,574, 7,928,727, 8,014,580, 8,058,601, 8,058,625, 8,086,011, 8,090,179, 8,094,898, and 8,103,487 for more complete discussions of the pixon method and its application.
FIG. 2 illustrates a generic imaging system 200 with an imaging detector 210, and a pixon reconstruction unit 220. The reconstruction is based on a pixon method that uses a pixon map P, which interacts with a pixon reconstruction algorithm 230. The pixon method refers to a method that smoothes each point in object space (hereafter an “object point”) by assigning a shape or volume to each object point as a basis for the pixon smoothing. The object space is the space in which the result of the image reconstruction is defined and which corresponds to a domain that was imaged using the imaging system 200. (It should be noted that “image space” is a synonymous term to “object space”, and the two terms are hereafter used interchangeably.) A corresponding data space is given by the data points measured with the imaging detector 210.
The pixon method provides high quality reconstruction of an image object I in object space from a measured data set d in data space. As a spatially adaptive reconstruction method, the pixon method applies a data-motivated smoothing operation to every object point. In doing so, the pixon method uses the principal of minimum complexity when assigning to every object point a pixon kernel function, which is the basis for the smoothing operation. Within the pixon reconstruction unit 220, the pixon map P defines which of the pixon kernel functions is assigned to each of the object points.
In the imaging system 200, the imaging detector 210 detects and communicates the measured data set d to the pixon reconstruction unit 220. The pixon reconstruction unit 220 uses the specially-adapted pixon reconstruction algorithms 230 to reconstruct the acquired data set d into an image object I. In doing so, the pixon reconstruction algorithm 230 uses a system matrix H to describe the properties of the imaging system 200, and to estimate an iteratively improved image object by adjusting the data model, which is the basis for the image object I. The image object I is, for example, displayed on a display 240 using well-known rendering techniques.
For every object point, the pixon map P provides a pixon kernel function that is determined on the basis of a minimum complexity method. This pixon kernel function is used in a pixon smoothing operation applied in object space.
The pixon method can also enable compressed sensing in the form of super-resolution, utilizing the normegativity of the image and minimum complexity to reconstruct an image with finer pixels than those with which the data are obtained. This is not a violation of the sampling theorem due to Nyquist (1928) and Shanon (1949) because of the extra conditions of normegativity and minimum complexity (e.g., Puetter et al 2005). Spatial frequencies beyond the diffraction limit, which are truncated in the data, can similarly be reconstructed in the image.
FIG. 3 illustrates an exemplary process flow of the pixon method. Pixon smoothing is applied sequentially to a standard reconstruction algorithm.
Using a standard reconstruction algorithm, the input image is fitted to a measured data set d (step 300). In accordance with the above discussed use of the pixon kernel operator K, the resulting estimate of the image is called a pseudoimage. The pixon map P is determined using the pseudoimage and the measured data set d (step 310). The pseudoimage is also the initial object for the pixon smoothing operation (step 320). During the pixon smoothing operation (step 320), each object point of the pseudoimage is smoothed over a pixon kernel function. (In some variations of existing pixon methods, the pixon map can also be updated in each iteration by computing from the updated image.
Iterative image reconstruction methods iteratively fit image models to measured data and thus minimize the effect of noise on the final image. The result of a reconstruction algorithm is an approximated image that is fit to the measured data set d according to the rules of the algorithm.
In the pixon method, an approximated image can be used as an input object for pixon smoothing, for pixon reconstruction, and for the determination of the pixon map.
The pixon method includes a search for the broadest possible pixon kernel functions at each point in object space that together support an adequate fit of an object to the measured data set d. In particular, the pixon map assigns to each object point a specific pixon kernel function.
There can be drawbacks to the step of first computing a pseudoimage from which the pixon map can then be determined. This process requires more computations and runs the risk of introducing artifacts into the pseudoimage, which can bias the pixon map and hence the final image reconstruction. Furthermore, the determination of the pixon map works less well if the transformation from object space to data space is nonlocal (Bhatnagar & Cornwell 2004). For example, in interferometry and magnetic-resonance imaging, the data are the Fourier transforms of the image (plus noise), with each Fourier wave (basis function of the image) spreading over the entire image. Another example is the partially nonlocal transformation in tomography. The data are 2D projections of a 3D image (plus noise), a transformation which is local transverse to the projection direction but nonlocal along the projection direction.
In view of the foregoing, there is a need for an improved method for determining the pixon map within the pixon method.